This regroups most of the figures in the thesis (appendix figures are not included, except for the ROC resample figures, the BAU simulated land use change and the radar graph for change in flow for both chapters).

Data prep

Cleaning and formatting

We start by preparing the data from the outputs of scripts 2.6, 7.1 and 8.1. The results table shown as output summarised all simulation efforts for the thesis.

# Land use data 
trans_df <- read_csv("../outputs/final/land_use_trans_df.csv")
vals_df <- read_csv("../outputs/final/land_use_vals_df.csv")

# Results dir + mun and mrc shapefile
mun <- st_read("../data/mun/munic_SHP_clean.shp", quiet = TRUE)
mrc <- st_read("../data/raw/vector/mrc_SHP/mrc_polygone.shp",  quiet = TRUE)
mrc.mont <- mrc %>% filter(MRS_NM_REG=="Montérégie") 
mrc.mont.reproj <- st_transform(mrc.mont, raster::crs(mun))

# Classes
classes <- read_csv("../config/rcl_tables/land_use/recode_table.csv")
forest_classes <- read_csv("../config/rcl_tables/land_use/recode_table_forest.csv") %>% 
  rename(new_code = "ID", new_class = "Name")

classes_unique <- unique(classes[, c("new_code", "new_class")]) %>% 
  bind_rows(forest_classes)

# results key
sce_code_vec <- 
  as.vector(t(outer(c("BAU-", "R-", "CorPr-", "R-CorPr-", "R(T)-CorPr-"), 
                    c("Hist", "Base", "RCP8"), paste0)))
sce_code_vec_run <- 
  rep(c("BAU", "R", "CorPr", "R-CorPr", "R(T)-CorPr"), each=3)

results_clean <- read_csv("../outputs/final/stsim_run_results.csv") %>% 
  dplyr::select(scenarioId, name) %>% 
  mutate(name =  gsub(.$name, pattern = " \\(.+\\)", replacement = "")) %>% 
  mutate(sce =  paste0("sce_", .$scenarioId)) %>% 
  mutate(chapter = c("none", rep("both", 3), rep("chap_1", 3), rep("chap_2", 9))) %>% 
  mutate(splitted = str_split(name, " \\| ")) %>% 
  mutate(climate = unlist(map(splitted, ~unlist(.x[2])))) %>% 
  mutate(run = unlist(map(splitted, ~unlist(.x[1])))) %>% 
  dplyr::select(-splitted) %>% 
  replace_na(list(climate = "none")) %>% 
  mutate(code = c("Control", sce_code_vec)) %>% 
  mutate(code_run = c("control", sce_code_vec_run))

# Final extracted datasets
df_final <- read_csv("../outputs/final/final_df_current_density_part1.csv") %>%
  bind_rows(read_csv("../outputs/final/final_df_current_density_part2.csv")) %>% 
  mutate(timestep = (timestep*10)+1990, source = "model", zone = as.character(zone))
df_final_origin <- read_csv("../outputs/final/final_df_origin_current_density.csv") %>%
  mutate(timestep = timestep*10+1990, source = "model")

# Stratum key
stratum_key <- read_csv("../config/stsim/SecondaryStratum.csv") %>%
  mutate(ID = as.factor(ID)) %>%
  rename(zone=ID, MUS_NM_MUN=Name) %>%
  mutate(zone = as.character(zone)) %>% 
  left_join(df_final, by="zone") %>%
  filter(MUS_NM_MUN!="Not Monteregie")

# Sum all municipalities
df_summarised <- df_final %>%
  group_by(sce, timestep, species, iteration) %>% 
  summarise(sum_cur = sum(current)) %>% ungroup %>%
  mutate(source = "model")
df_origin_summarised <- df_final_origin %>%
  group_by(timestep, species) %>% 
  summarise(sum_cur = sum(mean)) %>% ungroup %>% 
  mutate(sce = "sce_0", source = "observation")

# Joined dataset
joined <- full_join(df_summarised, df_origin_summarised, 
                    by=c("sce", "source", "species", "timestep", "sum_cur")) %>% 
  left_join(results_clean, by = "sce")  %>% 
  replace_na(list(climate = "none", run = "historic run")) %>% 
  # Make factors
  mutate(sce = as.factor(sce), run = as.factor(run), sum_cur = 10*(sum_cur),
         climate = factor(climate, levels = c("none", "baseline",
                                              "RCP 8.5")),
         run = factor(run, levels = c("historic run", "BAU run", 
                                      "BAU run + corrs protection", "BAU run + ref", 
                                      "BAU run + corrs protection + ref",
                                      "BAU run + corrs protection + ref TARGETED")))

# Histogram data 
hist_original <- read_csv("../outputs/final/final_values_output_original.csv")
hist_true <- read_csv("../outputs/final/final_values_output_TRUE.csv") %>% 
  mutate(sce = "TRUE")
histograms <- read_csv("../outputs/final/final_values_output.csv") %>% 
  bind_rows(hist_original) %>% 
  bind_rows(hist_true) %>% 
  left_join(results_clean, by="sce") 

# SURF Data
surf <- read_csv("../surf/surf_output.csv") %>% 
  mutate(timestep = timestep*10+1990) %>% 
  left_join(results_clean, by=c("scenario"="scenarioId")) %>% 
  replace_na(list(climate = "none", run = "historic run")) %>% 
  # Make factors
  mutate(sce = as.factor(sce), run = as.factor(run),
         climate = factor(climate, levels = c("none", "historic",
                                              "baseline", "RCP 8.5")),
         run = factor(run, levels = c("historic run", "BAU run", 
                                      "BAU run + corrs protection", "BAU run + ref", 
                                      "BAU run + corrs protection + ref",
                                      "BAU run + corrs protection + ref TARGETED")))

# Bar plot data
bar_data <- read_csv("../outputs/final/final_bar_plot_data.csv")

results_clean

Figures

Mean current flow change

Historic

Change in mean flow (in % of the the 1990 flow) between 1990 and 2010, contrasting observed historic change (1990-2010) and simulated change from historic model run.

the_width = 18 
the_height = 15
the_dpi = 500

fig_1_historic <- joined %>% 
  filter(climate == "none") %>% 
  group_by(scenarioId, species, iteration) %>% 
  group_modify(~mutate(.x , baseline = subset(.x, timestep == 1990)$sum_cur, 
                       the_diff = ((sum_cur-baseline)/baseline)*100 )) %>% 
  ggplot(aes(x=timestep, y=the_diff, col=source)) +
  scale_color_manual(values=c('#d8b365','#5ab4ac'), 
                     labels = c("Model", "Observation")) +
  geom_smooth(aes(group = sce), method = "glm", se=FALSE) +
  scale_linetype_manual(values = c(1,3,2,5))+
  geom_point(aes(group = sce)) +
  facet_wrap(~species, scales = "fixed") +
  labs(y = "Current flow (% Change)",
       x = "Year",
       col = "Source") +
  NULL
#fig_1_historic
ggsave(fig_1_historic,
       filename = "../thesis/figures/connectivity_decrease_x5species_historic.png",
       width = the_width, height = the_height, dpi = the_dpi)

Chap 1

Line graph

Change in mean flow (in % of the the 2010 flow) between 2010 (observed) and 2100 (predicted), contrasting BAU scenario (solid line) with other land use change scenarios - linear graph.

colors_sce <- 
  data.frame(sce = c("Baseline", "RCP85"), 
             color = c("#8DA0CB", "#800026"), stringsAsFactors = F)

fig_1_static_chap_1 <- joined %>% 
  filter(climate != "none", chapter %in% c("none", "both", "chap_1")) %>% 
  group_by(scenarioId, species, iteration) %>% 
  group_modify(~mutate(.x , baseline = subset(.x, timestep == 2010)$sum_cur, 
                       the_diff = ((sum_cur-baseline)/baseline)*100 )) %>% 
  ggplot(aes(x=timestep, y=the_diff, col=climate)) +
  scale_color_manual(values = colors_sce$color, 
                     labels = colors_sce$sce) +
  geom_smooth(aes(group = sce, linetype=code_run), alpha=0.5, method="loess") +
  #geom_point(aes(group = sce)) +
  scale_linetype_manual(values=c("solid", "twodash")) +
  facet_wrap(~species, scales = "fixed") +
  geom_hline(yintercept = 0, linetype = "dashed", alpha = 0.5) +
  #geom_vline(xintercept = 2020, linetype = "dashed", color = "darkGreen", alpha = 0.5) +
  labs(y = "Current flow (% Change)",
       x = "Year",
       col = "Climate",
       linetype = "Run") +
  guides(col = guide_legend(ncol = 2, nrow=2), 
         linetype = guide_legend(nrow=1, override.aes = list(colour = 'black'))) +
  theme(legend.position = c(0.85, 0.15))+
  NULL
#fig_1_static_chap_1
ggsave(fig_1_static_chap_1, 
       filename = "../thesis/figures/connectivity_decrease_x5species_chap1.png", 
       width = the_width, height = the_height, dpi = the_dpi)

Radar Graph

Change in mean flow (in % of the the 2010 flow) between 2010 (observed) and 2100 (predicted), contrasting BAU scenario with other land use change scenarios - radar graph.

# Default size
gridline.label.offset   = 4
legend.text.size = 20
group.line.width = 0.9
grid.label.size = 10
background.circle.colour    = "#ffffff"
group.point.size    = 5
axis.label.size = 8
group.colours = RColorBrewer::brewer.pal(name = "Paired",n=10)[c(2,4,6,8,10)]


radar_data_chap1 <- joined %>% 
  filter(climate != "none", chapter %in% c("none", "both", "chap_1")) %>% 
  group_by(timestep, species, code) %>% 
  summarise(sum_cur = mean(sum_cur)) %>% ungroup %>% 
  pivot_wider(names_from = timestep, values_from = sum_cur) %>% 
  mutate(diff = ((`2100`-`2010`)/`2010`)*100) %>% 
  select(-c(`2010`:`2100`)) %>%  
  pivot_wider(names_from = code, values_from = diff) %>% 
  rename(group = species)

radar_chap1 <-  
  ggradar(radar_data_chap1, centre.y = -20, legend.position = "right",
          grid.min = -20, grid.max = 3, grid.mid = 0, 
          values.radar = c("-20 %", "0 %", ""), 
          group.colours = group.colours,
          gridline.label.offset = gridline.label.offset, 
          legend.text.size = legend.text.size,
          group.line.width = group.line.width, 
          grid.label.size = grid.label.size,
          background.circle.colour = background.circle.colour, 
          group.point.size  = group.point.size, 
          axis.label.size = axis.label.size)
#radar_chap1
ggsave("../thesis/figures/radar_ggradar_chap1.png", radar_chap1,
       width = the_width, height = the_height, dpi = the_dpi)

Chap 2

Change in mean flow (in % of the the 2010 flow) between 2010 (observed) and 2100 (predicted), contrasting BAU scenario (solid line) with other conservation scenarios - linear graph.

fig_1_static_chap_2 <- joined %>% 
  filter(climate != "none", chapter %in% c("none", "both", "chap_2")) %>% 
  group_by(scenarioId, species, iteration) %>% 
  group_modify(~mutate(.x , baseline = subset(.x, timestep == 2010)$sum_cur, 
                       the_diff = ((sum_cur-baseline)/baseline)*100 )) %>% 
  ggplot(aes(x=timestep, y=the_diff, col=climate)) +
  scale_color_manual(values = colors_sce$color,
                     labels = colors_sce$sce) +
  geom_smooth(aes(group = sce, linetype=code_run), alpha=0.5, method="loess") +
  scale_linetype_manual(values=c("solid", "dashed", "longdash", "dotted")) +
  facet_wrap(~species, scales = "fixed") +
  geom_hline(yintercept = 0, linetype = "dashed", alpha = 0.5) +
  geom_vline(xintercept = 2020, linetype = "dashed", color = "darkGreen", alpha = 0.5) +
  labs(y = "Current flow (% Change)",
       x = "Year",
       col = "Climate",
       linetype = "Run") +
  guides(col = guide_legend(ncol = 2, nrow=2), 
         linetype = guide_legend(nrow = 2, ncol=2, 
                                 override.aes = list(colour = 'black'))) +
  theme(legend.position = c(0.85, 0.15))+
  NULL
#fig_1_static_chap_2
ggsave(fig_1_static_chap_2, 
       filename = "../thesis/figures/connectivity_decrease_x5species_chap2.png", 
       width = the_width, height = the_height, dpi = the_dpi)

Radar Graph

Change in mean flow (in % of the the 2010 flow) between 2010 (observed) and 2100 (predicted), contrasting BAU scenario (solid line) with other conservation scenarios - radar graph.

radar_data_chap2 <- joined %>% 
  filter(climate != "none", chapter %in% c("none", "both","chap_2")) %>% 
  group_by(timestep, species, code) %>% 
  summarise(sum_cur = mean(sum_cur)) %>% ungroup %>% 
  pivot_wider(names_from = timestep, values_from = sum_cur) %>% 
  mutate(diff = ((`2100`-`2010`)/`2010`)*100) %>% 
  select(-c(`2010`:`2100`)) %>%  
  pivot_wider(names_from = code, values_from = diff) %>% 
  rename(group = species)

radar_chap2 <-  
  ggradar(radar_data_chap2, centre.y = -20, legend.position = "right",
          grid.min = -20, grid.max = 3, grid.mid = 0, 
          values.radar = c("-20 %", "0 %", ""), 
          group.colours = group.colours,
          gridline.label.offset = gridline.label.offset, 
          legend.text.size = legend.text.size,
          group.line.width = group.line.width, 
          grid.label.size = grid.label.size,
          background.circle.colour = background.circle.colour, 
          group.point.size  = group.point.size, 
          axis.label.size = axis.label.size-2)
#radar_chap2
ggsave("../thesis/figures/radar_ggradar_chap2.png", radar_chap2,
       width = the_width, height = the_height, dpi = the_dpi)

Both

Radar Graph

Change in mean flow (in % of the the 2010 flow) between 2010 (observed) and 2100 (predicted), contrasting BAU scenario (solid line) with all other scenarios from borth chapters - radar graph.

radar_data_all <- joined %>% 
  filter(climate != "none") %>% 
  group_by(timestep, species, code) %>% 
  summarise(sum_cur = mean(sum_cur)) %>% ungroup %>% 
  pivot_wider(names_from = timestep, values_from = sum_cur) %>% 
  mutate(diff = ((`2100`-`2010`)/`2010`)*100) %>% 
  select(-c(`2010`:`2100`)) %>%  
  pivot_wider(names_from = code, values_from = diff) %>% 
  rename(group = species)

radar_all <-  
  ggradar(radar_data_all, centre.y = -20, legend.position = "right",
          grid.min = -20, grid.max = 3, grid.mid = 0, 
          values.radar = c("-20 %", "0 %", ""), 
          group.colours = group.colours,
          gridline.label.offset = gridline.label.offset, 
          legend.text.size = legend.text.size,
          group.line.width = group.line.width, 
          grid.label.size = grid.label.size,
          background.circle.colour = background.circle.colour, 
          group.point.size  = group.point.size, 
          axis.label.size = axis.label.size-2)
#radar_all
ggsave("../thesis/figures/radar_ggradar_both.png", radar_all,
       width = the_width, height = the_height, dpi = the_dpi)

ROC Curves

ROC curves from fitting both models (with resampling), plotted with the patchwork package.

# Roc curves data
urb <- readRDS("../data/temp/fit_rs_outcome_rf_urb_2.RDS")
agex <- readRDS("../data/temp/fit_rs_outcome_rf_agex_2.RDS")

p1 <- bind_rows(urb, .id = "fold") %>% 
  mutate(Fold = factor(fold, levels = as.character(1:10))) %>% 
  group_by(Fold) %>% 
  roc_curve(.pred, truth = outcome_fact) %>% 
  autoplot(add=T) +
  ggtitle(paste("Urbanisation")) +
  annotate(x = 0.75, y = 0.25, geom="label", size = 3,
           label = as.character(paste("Av AUC = 0.938 +/- 0.002"))) +
  theme(legend.position = "none") +
  NULL

p2 <- bind_rows(agex, .id = "fold") %>% 
  mutate(Fold = factor(fold, levels = as.character(1:10))) %>% 
  group_by(Fold) %>% 
  roc_curve(.pred, truth = outcome_fact) %>% 
  autoplot(add=T) +
  ggtitle(paste("Agricultural Expansion")) +
  annotate(x = 0.75, y = 0.25, geom="label", size = 3,
           label = as.character(paste("Av AUC = 0.929 +/- 0.002"))) +
  NULL

full = p1 + p2

#full
ggsave("../thesis/figures/double_roc_resample.png", full)

Histograms

Historic

Histograms of flow values for the year 2010 comparing the historic model run (1990-2010) with observations.

bi_colors <- c("#67a9cf", "#ef8a62")

hist_plot_origin <- histograms %>%
  filter(sce %in% c("TRUE", "sce_37")) %>% 
  filter(ts == 2) %>% 
  mutate(ts = ts*10+1990) %>% 
  mutate(ts = as.factor(ts)) %>% 
  ggplot(aes(x=bins, y=n, fill=sce, colour =sce)) + 
  geom_density(stat='identity', show.legend=T, alpha=0.3)+
  scale_fill_manual(values = bi_colors, name="Source", labels=c("Model", "Observation")) +
  scale_color_manual(values = bi_colors, name="Source", labels=c("Model", "Observation")) +
  facet_wrap(~species) +
  labs(x = "Flow intensity distribution (log)", 
       y = "") +
  theme(strip.text.y = element_text(angle=360, size=10, hjust = 0), 
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(), 
        axis.text.y.left = element_blank(),
        axis.ticks.y.left = element_blank()) +
  NULL
#hist_plot_origin
ggsave("../thesis/figures/original_hists.png", hist_plot_origin,
       width = the_width, height = the_height, dpi = the_dpi)

Chap 1

Histograms of flow values for all species and scenarios, for 2010 (observed) and 2100 (predicted).

bi_colors <- c("#67a9cf", "#ef8a62")

hist_plot_1 <- histograms %>% 
  mutate(ts = as.factor(ts)) %>% 
  filter(chapter %in% c('both', 'chap_1')) %>%  
  ggplot(aes(x=bins, y=n, group=ts)) + 
  geom_density(stat='identity', show.legend=T, aes(fill=factor(ts), color=factor(ts)), alpha=0.3)+
  facet_grid(code~species) +
  scale_fill_manual(values = bi_colors, 
                    labels = c("2010", "2100")) +
  scale_color_manual(values = bi_colors, 
                     labels = c("2010", "2100")) +
  labs(x = "Flow intensity distribution (log)", 
       y = "")+
  theme(strip.text.y = element_text(angle=360, size=10, hjust = 0), 
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(), 
        axis.text.y.left = element_blank(),
        axis.ticks.y.left = element_blank(), 
        legend.position = c(0.95, 0.96)) +
  labs(color="Timestep", fill="Timestep") +
  NULL
#hist_plot_1
ggsave("../thesis/figures/hist_chap1.png", hist_plot_1, width = 18, height = 15)

Chap 2

Histograms of flow values for all species and scenarios, for 2010 (observed) and 2100 (predicted).

hist_plot_2 <- histograms %>% 
  mutate(ts = as.factor(ts)) %>% 
  filter(chapter %in% c('chap_2')) %>%  
  ggplot(aes(x=bins, y=n, group=ts)) + 
  geom_density(stat='identity', aes(fill=factor(ts),  colour=factor(ts)), alpha=0.3)+
  facet_grid(code~species) +
  scale_fill_manual(values = bi_colors, 
                    labels = c("2010", "2100")) +
  scale_color_manual(values = bi_colors, 
                     labels = c("2010", "2100")) +
  labs(x = "Flow intensity distribution (log)", 
       y = "") +
  theme(strip.text.y = element_text(angle=360, size=10, hjust = 0), 
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(), 
        axis.text.y.left = element_blank(),
        axis.ticks.y.left = element_blank(),
        legend.position = c(0.95, 0.96)) +
  labs(color="Timestep", fill="Timestep") +
  NULL
#hist_plot_2
ggsave("../thesis/figures/hist_chap2.png", hist_plot_2, width = 18, height = 20)

SURF Analysis

Chap1

Linear Graph

Change in the number of features (in % of the the 2010 flow) identified by the SURF analyssis between 2010 (observed) and 2100 (predicted), contrasting BAU scenario with other land use change scenarios - linear graph.

surf_1 <- surf %>% 
  mutate(scenario = as.factor(scenario)) %>% 
  filter(climate != "none", chapter %in% c("none", "both", "chap_1")) %>% 
  filter(climate != "none") %>% 
  group_by(scenario, species, iter) %>% 
  group_modify(~mutate(.x , baseline = subset(.x, timestep == 2010)$kp_nb, 
                       the_diff = ((kp_nb-baseline)/baseline)*100 )) %>% ungroup %>% 
  ggplot(aes(x = timestep, y = the_diff, color=climate)) +
  scale_linetype_manual(values=c("solid", "twodash")) +
  geom_smooth(aes(group = scenario, linetype=code_run), alpha=0.2, method="loess")+
  #geom_point(aes(group=scenario), alpha= 0.05, size = 0.5)+
  facet_wrap(~species, scales = "fixed")+
  # scale_color_manual(values = colors_sce$color,
  #                    labels = colors_sce$sce) +
  labs(y = "Change in nb of Keypoints detected (%)",
       x = "Year",
       col = "Climate",
       linetype = "Run") +
  guides(col = guide_legend(ncol = 2, nrow=2), 
         linetype = guide_legend(nrow=1, override.aes = list(colour = 'black'))) +
  NULL
#surf_1
ggsave("../thesis/figures/surf_chap1.png", surf_1,
       width = the_width, height = the_height, dpi = the_dpi)

Radar Graph

Change in the number of features (in % of the the 2010 flow) identified by the SURF analyssis between 2010 (observed) and 2100 (predicted), contrasting BAU scenario with other land use change scenarios - radar graph.

surf_radar_data_1 <- surf %>% 
  filter(climate != "none", chapter %in% c("none", "both","chap_1")) %>% 
  filter(sce != "sce_0", name != "historic run") %>% 
  group_by(timestep, species, code) %>% 
  summarise(kp_nb = mean(kp_nb)) %>% ungroup %>% 
  pivot_wider(names_from = timestep, values_from = kp_nb) %>% 
  mutate(diff = ((`2100`-`2010`)/`2010`)*100) %>% 
  select(-c(`2010`:`2100`)) %>%  
  pivot_wider(names_from = code, values_from = diff) %>% 
  rename(group = species)

surf_radar_1 <- 
  ggradar(surf_radar_data_1, centre.y = -60, legend.position = "right",
          grid.min = -60, grid.max = 5, grid.mid = 0, 
          values.radar = c("-60 %", "0 %", ""), 
          group.colours = group.colours,
          gridline.label.offset = 20, 
          legend.text.size = legend.text.size,
          group.line.width = group.line.width, 
          grid.label.size = grid.label.size,
          background.circle.colour = background.circle.colour, 
          group.point.size  = group.point.size, 
          axis.label.size = axis.label.size)
# surf_radar_1
ggsave("../thesis/figures/surf_radar_chap1.png", surf_radar_1,
       width = the_width, height = the_height, dpi = the_dpi)

Chap 2

Linear Graph

Change in the number of features (in % of the the 2010 flow) identified by the SURF analyssis between 2010 (observed) and 2100 (predicted), contrasting BAU scenario with other conservation scenarios - linear graph.

surf_2 <- surf %>% 
  mutate(scenario = as.factor(scenario)) %>% 
  filter(climate != "none", chapter %in% c("none", "both", "chap_2")) %>% 
  filter(climate != "none") %>% 
  group_by(scenario, species, iter) %>% 
  group_modify(~mutate(.x , baseline = subset(.x, timestep == 2010)$kp_nb, 
                       the_diff = ((kp_nb-baseline)/baseline)*100 )) %>% ungroup %>% 
  ggplot(aes(x = timestep, y = the_diff, color=climate)) +
  geom_smooth(aes(group = scenario, linetype=code_run), alpha=0.2, method="loess")+
  #geom_point(aes(group=scenario), alpha= 0.05, size = 0.5)+
  scale_linetype_manual(values=c("solid", "dashed", "longdash", "dotted")) +
  facet_wrap(~species, scales = "fixed")+
  # scale_color_manual(values = colors_sce$color,
  #                    labels = colors_sce$sce) +
  labs(y = "Change in nb of Keypoints detected (%)",
       x = "Year",
       col = "Climate",
       linetype = "Run") +
  guides(col = guide_legend(ncol = 2, nrow=2), 
         linetype = guide_legend(nrow=2, ncol=2, 
                                 override.aes = list(colour = 'black'))) +
  NULL
#surf_2
ggsave("../thesis/figures/surf_chap2.png", surf_2,
       width = the_width, height = the_height, dpi = the_dpi)

Radar Graph

Change in the number of features (in % of the the 2010 flow) identified by the SURF analyssis between 2010 (observed) and 2100 (predicted), contrasting BAU scenario with other conservation scenarios - radar graph.

surf_radar_data_2 <- surf %>% 
  filter(climate != "none", chapter %in% c("none", "both","chap_2")) %>% 
  filter(sce != "sce_0", name != "historic run") %>% 
  group_by(timestep, species, code) %>% 
  summarise(kp_nb = mean(kp_nb)) %>% ungroup %>% 
  pivot_wider(names_from = timestep, values_from = kp_nb) %>% 
  mutate(diff = ((`2100`-`2010`)/`2010`)*100) %>% 
  select(-c(`2010`:`2100`)) %>%  
  pivot_wider(names_from = code, values_from = diff) %>% 
  rename(group = species)

surf_radar_2 <- 
  ggradar(surf_radar_data_2, centre.y = -60, legend.position = "right",
          grid.min = -60, grid.max = 5, grid.mid = 0, 
          values.radar = c("-60 %", "0 %", ""), 
          group.colours = group.colours,
          gridline.label.offset = 22, 
          legend.text.size = legend.text.size,
          group.line.width = group.line.width, 
          grid.label.size = grid.label.size,
          background.circle.colour = background.circle.colour, 
          group.point.size  = group.point.size, 
          axis.label.size = axis.label.size-2)
#surf_radar_2
ggsave("../thesis/figures/surf_radar_chap2.png", surf_radar_2, 
       width = the_width, height = the_height, dpi = the_dpi)

Prediction maps

list_of_files <-  list.files("../libraries/stsim/monteregie-conncons-scripted.ssim.output/Scenario-37/stsim_OutputSpatialState", full.names = T)
list_of_files <- list_of_files[grep(x=list_of_files, "ts2")]
stratum <- raster("../data/stsim/aggregated/primary_stratum_mont_or_not_or_PA.tif")

Predicted probability of urbanization (left) vs observed urbanization (right) between 1990 and 2010.

names <- paste0(c("historic_compare_agex", "historic_compare_urb"), ".png")

for (class in c(1:2)){
  
  the_stack <- stack(lapply(list_of_files, raster)) == class
  original <- raster("../data/land_use/aggregated/aggregated_lu_buffered_1990_patched.tif") == class
  final <- raster("../data/land_use/aggregated/aggregated_lu_buffered_2010_patched.tif") == class
  
  # Transform, give names
  the_stack[original == 1] <- 0
  final[original == 1] <- 0
  the_stack[stratum==0] <- NA
  final[stratum==0] <- NA
  the_stack <- trim(the_stack)
  final <- trim(final)
  
  the_mean <-  mean(the_stack)
  names(the_mean) <- "mean"
  names(final) <- "observation"
  
  # test <- the_mean
  # test[(which(values(!is.na(final)) & values(is.na(the_mean))))] <- 100
  # plot(test)
  
  # Fill the "other" values
  the_mean[(which(values(!is.na(final)) & values(is.na(the_mean))))] <- 0
  
  stack_to_plot <- stack(the_mean, final) 
  names(stack_to_plot) <- c("Model", "Observations")
  
  mxp <- 1e10000 # 1000
  
  myPal <- RColorBrewer::brewer.pal('OrRd', n=9)
  myTheme <- rasterTheme(region = myPal)
  probs <- levelplot(stack_to_plot$Model, scales=list(draw=FALSE),
                     maxpixels = mxp, 
                     par.settings = myTheme,
                     margin=FALSE, auto.key=FALSE, 
                     colorkey=list(space="bottom"), 
                     main = "A)")
  probs <- as.ggplot(probs) 
  
  stack_to_plot$Observations <- ratify(stack_to_plot$Observations)
  rat <- levels(stack_to_plot$Observations)[[1]]
  rat$Urbanised <- c("Not Observed", "Observed")
  levels(stack_to_plot$Observations) <- rat
  bin <- levelplot(stack_to_plot$Observations, col.regions=c('grey', 'indianred4'), 
                   scales=list(draw=FALSE), 
                   maxpixels = mxp, 
                   colorkey=list(space="bottom", height=0.3), 
                   main = "B)")
  bin <- as.ggplot(bin)
  
  final_plot <- probs + bin
  
  ggsave(plot = final_plot, filename = paste0("../thesis/figures/",names[class]), width = 20, height = 10, dpi = the_dpi)
}

Predicted probability of agricultural expansion (left) vs observed agricultural expansion (right) between 1990 and 2010.

Land use maps

Land use maps showing predictions of land use changes (Only BAU is displayed below).

source('../scripts/functions/draw_scenario.R')
for (scenario in c(39, 42, 45, 48, 51)){
  the_plot <- draw_scenario(scenario, mxp = mxp)
  the_name <- paste0("scenario_", as.character(scenario),"_compare.png")
  ggsave(plot = the_plot, 
         filename = paste0("../thesis/figures/",the_name), width = 30, height = 12, dpi = the_dpi)
}

Comparison of Monteregie at the beginning of the BAU scenario in 2010, and at the end in 2100 (under baseline climate scenario).

Bar charts

Both scenarios

Land use class frequencies in 2010 (observed) and 2100 (predicted) for the BAU and Reforestation scenarios (under baseline climate scenario).

final_df_39 <- final_df_39 %>% 
  mutate(sce = "BAU") %>% 
  mutate(type = ifelse(value <10, "non_forest", "forest"))
final_df_42 <- final_df_42 %>% 
  mutate(sce = "Reforestation")  %>% 
  mutate(type = ifelse(value <10, "non_forest", "forest"))

final_all <- bind_rows(final_df_39, final_df_42)

forest <- final_all %>% 
  filter(type == "forest") %>% 
  mutate(value = as.factor(value), timestep = as.factor(timestep)) %>% 
  ggplot(aes(fill=timestep, y=area, x=new_class)) + 
  geom_bar(position="dodge", stat="identity") +
  scale_fill_manual(values = bi_colors) +
  theme(legend.position = "none") +
  facet_grid(~sce) +
  labs(x = "Land use Class", 
       y = "Mean Area (hectare)", 
       fill = "Year") +
  NULL

non_forest <- final_all %>% 
  filter(type != "forest") %>% 
  filter(!(new_class %in% c("Water", "Wetlands", "Roads"))) %>% 
  mutate(value = as.factor(value), timestep = as.factor(timestep)) %>% 
  ggplot(aes(fill=timestep, y=area, x=new_class)) + 
  geom_bar(position="dodge", stat="identity") +
  scale_fill_manual(values = bi_colors) +
  theme(legend.position = "right") +
  facet_grid(~sce) +
  labs(x = "Land use Class", 
       y = "Mean Area (hectare)", 
       fill = "Year") +
  NULL

assembled <- wrap_plots(forest, non_forest, ncol=1, nrow=2)
#assembled
ggsave("../thesis/figures/bar_both.png", assembled,
       width = the_width, height = 20, dpi = the_dpi)

Land use class frequencies in 2010 (observed) and 2100 (predicted) for the BAU and conservation scenarios (under baseline climate scenario).

Networks

CHAP 1
networks <- read_csv("../outputs/final/final_df_network_stats.csv") %>% 
  left_join(results_clean, by = c("scenario" = "scenarioId"))

networks %>% 
  mutate(timestep = timestep*10 + 1990) %>% 
  filter(climate != "none", chapter %in% c("none", "both", "chap_1")) %>% 
  group_by(scenario, species, iteration) %>% 
  group_modify(~mutate(.x , baseline = subset(.x, timestep == 2010)$overallECgap, 
                       the_diff = ((overallECgap-baseline)/baseline)*100 )) %>%
  ggplot(aes(x=timestep, y=the_diff, col=climate)) +
  scale_color_manual(values = colors_sce$color, 
                     labels = colors_sce$sce) +
  geom_smooth(aes(group = sce, linetype=code_run), alpha=0.5, method="loess") +
  #geom_point(aes(group = sce)) +
  scale_linetype_manual(values=c("solid", "twodash")) +
  facet_wrap(~species, scales = "fixed") +
  geom_hline(yintercept = 0, linetype = "dashed", alpha = 0.5) +
  #geom_vline(xintercept = 2020, linetype = "dashed", color = "darkGreen", alpha = 0.5) +
  labs(y = "overallECgap",
       x = "Year",
       col = "Climate",
       linetype = "Run") +
  guides(col = guide_legend(ncol = 2, nrow=2), 
         linetype = guide_legend(nrow=1, override.aes = list(colour = 'black'))) +
  theme(legend.position = c(0.85, 0.15))+
  NULL -> net_1

ggsave(net_1, 
       filename = "../thesis/figures/ECGAP.png", 
       width = the_width, height = the_height, dpi = the_dpi)
networks %>% 
  mutate(timestep = timestep*10 + 1990) %>% 
  filter(climate != "none", chapter %in% c("none", "both", "chap_1")) %>% 
  group_by(scenario, species, iteration) %>% 
  group_modify(~mutate(.x , baseline = subset(.x, timestep == 2010)$overallECnatal, 
                       the_diff = ((overallECnatal-baseline)/baseline)*100 )) %>%
  ggplot(aes(x=timestep, y=the_diff, col=climate)) +
  scale_color_manual(values = colors_sce$color, 
                     labels = colors_sce$sce) +
  geom_smooth(aes(group = sce, linetype=code_run), alpha=0.5, method="loess") +
  #geom_point(aes(group = sce)) +
  scale_linetype_manual(values=c("solid", "twodash")) +
  facet_wrap(~species, scales = "fixed") +
  geom_hline(yintercept = 0, linetype = "dashed", alpha = 0.5) +
  #geom_vline(xintercept = 2020, linetype = "dashed", color = "darkGreen", alpha = 0.5) +
  labs(y = "overallECnatal",
       x = "Year",
       col = "Climate",
       linetype = "Run") +
  guides(col = guide_legend(ncol = 2, nrow=2), 
         linetype = guide_legend(nrow=1, override.aes = list(colour = 'black'))) +
  theme(legend.position = c(0.85, 0.15))+
  NULL -> net_2

ggsave(net_2, 
       filename = "../thesis/figures/ECNATAL.png", 
       width = the_width, height = the_height, dpi = the_dpi)
networks %>% 
  mutate(timestep = timestep*10 + 1990) %>% 
  filter(climate != "none", chapter %in% c("none", "both", "chap_1")) %>% 
  group_by(scenario, species, iteration) %>% 
  group_modify(~mutate(.x , baseline = subset(.x, timestep == 2010)$numNodes, 
                       the_diff = ((numNodes-baseline)/baseline)*100 )) %>%
  ggplot(aes(x=timestep, y=the_diff, col=climate)) +
  scale_color_manual(values = colors_sce$color, 
                     labels = colors_sce$sce) +
  geom_smooth(aes(group = sce, linetype=code_run), alpha=0.5, method="loess") +
  #geom_point(aes(group = sce)) +
  scale_linetype_manual(values=c("solid", "twodash")) +
  facet_wrap(~species, scales = "fixed") +
  geom_hline(yintercept = 0, linetype = "dashed", alpha = 0.5) +
  #geom_vline(xintercept = 2020, linetype = "dashed", color = "darkGreen", alpha = 0.5) +
  labs(y = "numNodes",
       x = "Year",
       col = "Climate",
       linetype = "Run") +
  guides(col = guide_legend(ncol = 2, nrow=2), 
         linetype = guide_legend(nrow=1, override.aes = list(colour = 'black'))) +
  theme(legend.position = c(0.85, 0.15))+
  NULL -> net_3

ggsave(net_3, 
       filename = "../thesis/figures/NUMNODES.png", 
       width = the_width, height = the_height, dpi = the_dpi)
networks %>% 
  mutate(timestep = timestep*10 + 1990) %>% 
  filter(climate != "none", chapter %in% c("none", "both", "chap_1")) %>% 
  group_by(scenario, species, iteration) %>% 
  group_modify(~mutate(.x , baseline = subset(.x, timestep == 2010)$numLinks, 
                       the_diff = ((numLinks-baseline)/baseline)*100 )) %>%
  ggplot(aes(x=timestep, y=the_diff, col=climate)) +
  scale_color_manual(values = colors_sce$color, 
                     labels = colors_sce$sce) +
  geom_smooth(aes(group = sce, linetype=code_run), alpha=0.5, method="loess") +
  #geom_point(aes(group = sce)) +
  scale_linetype_manual(values=c("solid", "twodash")) +
  facet_wrap(~species, scales = "fixed") +
  geom_hline(yintercept = 0, linetype = "dashed", alpha = 0.5) +
  #geom_vline(xintercept = 2020, linetype = "dashed", color = "darkGreen", alpha = 0.5) +
  labs(y = "numLinks",
       x = "Year",
       col = "Climate",
       linetype = "Run") +
  guides(col = guide_legend(ncol = 2, nrow=2), 
         linetype = guide_legend(nrow=1, override.aes = list(colour = 'black'))) +
  theme(legend.position = c(0.85, 0.15))+
  NULL -> net_4

ggsave(net_4, 
       filename = "../thesis/figures/NUMLINKS.png", 
       width = the_width, height = the_height, dpi = the_dpi)
networks %>% 
  mutate(timestep = timestep*10 + 1990) %>% 
  filter(climate != "none", chapter %in% c("none", "both", "chap_1")) %>% 
  group_by(scenario, species, iteration) %>% 
  group_modify(~mutate(.x , baseline = subset(.x, timestep == 2010)$meanAreaQuality, 
                       the_diff = ((meanAreaQuality-baseline)/baseline)*100 )) %>%
  ggplot(aes(x=timestep, y=the_diff, col=climate)) +
  scale_color_manual(values = colors_sce$color, 
                     labels = colors_sce$sce) +
  geom_smooth(aes(group = sce, linetype=code_run), alpha=0.5, method="loess") +
  #geom_point(aes(group = sce)) +
  scale_linetype_manual(values=c("solid", "twodash")) +
  facet_wrap(~species, scales = "fixed") +
  geom_hline(yintercept = 0, linetype = "dashed", alpha = 0.5) +
  #geom_vline(xintercept = 2020, linetype = "dashed", color = "darkGreen", alpha = 0.5) +
  labs(y = "meanAreaQuality",
       x = "Year",
       col = "Climate",
       linetype = "Run") +
  guides(col = guide_legend(ncol = 2, nrow=2), 
         linetype = guide_legend(nrow=1, override.aes = list(colour = 'black'))) +
  theme(legend.position = c(0.85, 0.15))+
  NULL -> net_5

ggsave(net_5, 
       filename = "../thesis/figures/meanAreaQuality.png", 
       width = the_width, height = the_height, dpi = the_dpi)
networks %>% 
  mutate(timestep = timestep*10 + 1990) %>% 
  filter(climate != "none", chapter %in% c("none", "both", "chap_1")) %>% 
  group_by(scenario, species, iteration) %>% 
  group_modify(~mutate(.x , baseline = subset(.x, timestep == 2010)$meanCost, 
                       the_diff = ((meanCost-baseline)/baseline)*100 )) %>%
  ggplot(aes(x=timestep, y=the_diff, col=climate)) +
  scale_color_manual(values = colors_sce$color, 
                     labels = colors_sce$sce) +
  geom_smooth(aes(group = sce, linetype=code_run), alpha=0.5, method="loess") +
  #geom_point(aes(group = sce)) +
  scale_linetype_manual(values=c("solid", "twodash")) +
  facet_wrap(~species, scales = "fixed") +
  geom_hline(yintercept = 0, linetype = "dashed", alpha = 0.5) +
  #geom_vline(xintercept = 2020, linetype = "dashed", color = "darkGreen", alpha = 0.5) +
  labs(y = "meanCost",
       x = "Year",
       col = "Climate",
       linetype = "Run") +
  guides(col = guide_legend(ncol = 2, nrow=2), 
         linetype = guide_legend(nrow=1, override.aes = list(colour = 'black'))) +
  theme(legend.position = c(0.85, 0.15))+
  NULL -> net_6

ggsave(net_6, 
       filename = "../thesis/figures/meanCost.png", 
       width = the_width, height = the_height, dpi = the_dpi)
CHAP 2
networks %>% 
  mutate(timestep = timestep*10 + 1990) %>% 
  filter(climate != "none", chapter %in% c("none", "both", "chap_2")) %>% 
  group_by(scenario, species, iteration) %>% 
  group_modify(~mutate(.x , baseline = subset(.x, timestep == 2010)$overallECgap, 
                       the_diff = ((overallECgap-baseline)/baseline)*100 )) %>% 
  ggplot(aes(x=timestep, y=the_diff, col=climate)) +
  scale_color_manual(values = colors_sce$color,
                     labels = colors_sce$sce) +
  geom_smooth(aes(group = sce, linetype=code_run), alpha=0.5, method="loess") +
  scale_linetype_manual(values=c("solid", "dashed", "longdash", "dotted")) +
  facet_wrap(~species, scales = "fixed") +
  geom_hline(yintercept = 0, linetype = "dashed", alpha = 0.5) +
  geom_vline(xintercept = 2020, linetype = "dashed", color = "darkGreen", alpha = 0.5) +
  labs(y = "overallECgap",
       x = "Year",
       col = "Climate",
       linetype = "Run") +
  guides(col = guide_legend(ncol = 2, nrow=2), 
         linetype = guide_legend(nrow = 2, ncol=2, 
                                 override.aes = list(colour = 'black'))) +
  theme(legend.position = c(0.85, 0.15))+
  NULL -> net_1_2

ggsave(net_1_2, 
       filename = "../thesis/figures/ECGAP2.png", 
       width = the_width, height = the_height, dpi = the_dpi)
networks %>% 
  mutate(timestep = timestep*10 + 1990) %>% 
  filter(climate != "none", chapter %in% c("none", "both", "chap_2")) %>% 
  group_by(scenario, species, iteration) %>% 
  group_modify(~mutate(.x , baseline = subset(.x, timestep == 2010)$overallECnatal, 
                       the_diff = ((overallECnatal-baseline)/baseline)*100 )) %>% 
  ggplot(aes(x=timestep, y=the_diff, col=climate)) +
  scale_color_manual(values = colors_sce$color,
                     labels = colors_sce$sce) +
  geom_smooth(aes(group = sce, linetype=code_run), alpha=0.5, method="loess") +
  scale_linetype_manual(values=c("solid", "dashed", "longdash", "dotted")) +
  facet_wrap(~species, scales = "fixed") +
  geom_hline(yintercept = 0, linetype = "dashed", alpha = 0.5) +
  geom_vline(xintercept = 2020, linetype = "dashed", color = "darkGreen", alpha = 0.5) +
  labs(y = "overallECnatal",
       x = "Year",
       col = "Climate",
       linetype = "Run") +
  guides(col = guide_legend(ncol = 2, nrow=2), 
         linetype = guide_legend(nrow = 2, ncol=2, 
                                 override.aes = list(colour = 'black'))) +
  theme(legend.position = c(0.85, 0.15))+
  NULL -> net_2_2

ggsave(net_2_2, 
       filename = "../thesis/figures/ECNATAL2.png", 
       width = the_width, height = the_height, dpi = the_dpi)
networks %>% 
  mutate(timestep = timestep*10 + 1990) %>% 
  filter(climate != "none", chapter %in% c("none", "both", "chap_2")) %>% 
  group_by(scenario, species, iteration) %>% 
  group_modify(~mutate(.x , baseline = subset(.x, timestep == 2010)$numNodes, 
                       the_diff = ((numNodes-baseline)/baseline)*100 )) %>% 
  ggplot(aes(x=timestep, y=the_diff, col=climate)) +
  scale_color_manual(values = colors_sce$color,
                     labels = colors_sce$sce) +
  geom_smooth(aes(group = sce, linetype=code_run), alpha=0.5, method="loess") +
  scale_linetype_manual(values=c("solid", "dashed", "longdash", "dotted")) +
  facet_wrap(~species, scales = "fixed") +
  geom_hline(yintercept = 0, linetype = "dashed", alpha = 0.5) +
  geom_vline(xintercept = 2020, linetype = "dashed", color = "darkGreen", alpha = 0.5) +
  labs(y = "numNodes",
       x = "Year",
       col = "Climate",
       linetype = "Run") +
  guides(col = guide_legend(ncol = 2, nrow=2), 
         linetype = guide_legend(nrow = 2, ncol=2, 
                                 override.aes = list(colour = 'black'))) +
  theme(legend.position = c(0.85, 0.15))+
  NULL -> net_3_2

ggsave(net_3_2, 
       filename = "../thesis/figures/NUMNODES2.png", 
       width = the_width, height = the_height, dpi = the_dpi)
networks %>% 
  mutate(timestep = timestep*10 + 1990) %>% 
  filter(climate != "none", chapter %in% c("none", "both", "chap_2")) %>% 
  group_by(scenario, species, iteration) %>% 
  group_modify(~mutate(.x , baseline = subset(.x, timestep == 2010)$numLinks, 
                       the_diff = ((numLinks-baseline)/baseline)*100 )) %>% 
  ggplot(aes(x=timestep, y=the_diff, col=climate)) +
  scale_color_manual(values = colors_sce$color,
                     labels = colors_sce$sce) +
  geom_smooth(aes(group = sce, linetype=code_run), alpha=0.5, method="loess") +
  scale_linetype_manual(values=c("solid", "dashed", "longdash", "dotted")) +
  facet_wrap(~species, scales = "fixed") +
  geom_hline(yintercept = 0, linetype = "dashed", alpha = 0.5) +
  geom_vline(xintercept = 2020, linetype = "dashed", color = "darkGreen", alpha = 0.5) +
  labs(y = "numLinks",
       x = "Year",
       col = "Climate",
       linetype = "Run") +
  guides(col = guide_legend(ncol = 2, nrow=2), 
         linetype = guide_legend(nrow = 2, ncol=2, 
                                 override.aes = list(colour = 'black'))) +
  theme(legend.position = c(0.85, 0.15))+
  NULL -> net_4_2

ggsave(net_4_2, 
       filename = "../thesis/figures/NUMLINKS2.png", 
       width = the_width, height = the_height, dpi = the_dpi)
networks %>% 
  mutate(timestep = timestep*10 + 1990) %>% 
  filter(climate != "none", chapter %in% c("none", "both", "chap_2")) %>% 
  group_by(scenario, species, iteration) %>% 
  group_modify(~mutate(.x , baseline = subset(.x, timestep == 2010)$meanAreaQuality, 
                       the_diff = ((meanAreaQuality-baseline)/baseline)*100 )) %>% 
  ggplot(aes(x=timestep, y=the_diff, col=climate)) +
  scale_color_manual(values = colors_sce$color,
                     labels = colors_sce$sce) +
  geom_smooth(aes(group = sce, linetype=code_run), alpha=0.5, method="loess") +
  scale_linetype_manual(values=c("solid", "dashed", "longdash", "dotted")) +
  facet_wrap(~species, scales = "fixed") +
  geom_hline(yintercept = 0, linetype = "dashed", alpha = 0.5) +
  geom_vline(xintercept = 2020, linetype = "dashed", color = "darkGreen", alpha = 0.5) +
  labs(y = "meanAreaQuality",
       x = "Year",
       col = "Climate",
       linetype = "Run") +
  guides(col = guide_legend(ncol = 2, nrow=2), 
         linetype = guide_legend(nrow = 2, ncol=2, 
                                 override.aes = list(colour = 'black'))) +
  theme(legend.position = c(0.85, 0.15))+
  NULL -> net_5_2

ggsave(net_5_2, 
       filename = "../thesis/figures/meanAreaQuality_2.png", 
       width = the_width, height = the_height, dpi = the_dpi)
networks %>% 
  mutate(timestep = timestep*10 + 1990) %>% 
  filter(climate != "none", chapter %in% c("none", "both", "chap_2")) %>% 
  group_by(scenario, species, iteration) %>% 
  group_modify(~mutate(.x , baseline = subset(.x, timestep == 2010)$meanCost, 
                       the_diff = ((meanCost-baseline)/baseline)*100 )) %>% 
  ggplot(aes(x=timestep, y=the_diff, col=climate)) +
  scale_color_manual(values = colors_sce$color,
                     labels = colors_sce$sce) +
  geom_smooth(aes(group = sce, linetype=code_run), alpha=0.5, method="loess") +
  scale_linetype_manual(values=c("solid", "dashed", "longdash", "dotted")) +
  facet_wrap(~species, scales = "fixed") +
  geom_hline(yintercept = 0, linetype = "dashed", alpha = 0.5) +
  geom_vline(xintercept = 2020, linetype = "dashed", color = "darkGreen", alpha = 0.5) +
  labs(y = "meanCost",
       x = "Year",
       col = "Climate",
       linetype = "Run") +
  guides(col = guide_legend(ncol = 2, nrow=2), 
         linetype = guide_legend(nrow = 2, ncol=2, 
                                 override.aes = list(colour = 'black'))) +
  theme(legend.position = c(0.85, 0.15))+
  NULL -> net_6_2

ggsave(net_6_2, 
       filename = "../thesis/figures/meanCost_2.png", 
       width = the_width, height = the_height, dpi = the_dpi)